{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a05045a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rdata\n",
    "import os\n",
    "import fastkde\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4b1ede45",
   "metadata": {},
   "outputs": [
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'v11acs_res/v11acs_mag_pos.rda'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "Input \u001b[0;32mIn [2]\u001b[0m, in \u001b[0;36m<cell line: 3>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m nm_pos \u001b[38;5;241m=\u001b[39m rdata\u001b[38;5;241m.\u001b[39mparser\u001b[38;5;241m.\u001b[39mparse_file(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mv11acs_nm_pos.rda\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m      2\u001b[0m col_pos \u001b[38;5;241m=\u001b[39m rdata\u001b[38;5;241m.\u001b[39mparser\u001b[38;5;241m.\u001b[39mparse_file(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mv11acs_col_pos.rda\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m----> 3\u001b[0m mag_pos \u001b[38;5;241m=\u001b[39m \u001b[43mrdata\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mparser\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mparse_file\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mv11acs_res/v11acs_mag_pos.rda\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m      4\u001b[0m col_err_pos \u001b[38;5;241m=\u001b[39m rdata\u001b[38;5;241m.\u001b[39mparser\u001b[38;5;241m.\u001b[39mparse_file(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mv11acs_res/v11acs_col_err_pos.rda\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m      6\u001b[0m nm_pos \u001b[38;5;241m=\u001b[39m rdata\u001b[38;5;241m.\u001b[39mconversion\u001b[38;5;241m.\u001b[39mconvert(nm_pos)\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.9/site-packages/rdata/parser/_parser.py:1134\u001b[0m, in \u001b[0;36mparse_file\u001b[0;34m(file_or_path, expand_altrep, altrep_constructor_dict, extension)\u001b[0m\n\u001b[1;32m   1132\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m extension \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m   1133\u001b[0m         extension \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mgetattr\u001b[39m(path, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124msuffix\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28;01mNone\u001b[39;00m)\n\u001b[0;32m-> 1134\u001b[0m     data \u001b[38;5;241m=\u001b[39m \u001b[43mpath\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mread_bytes\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1136\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m parse_data(\n\u001b[1;32m   1137\u001b[0m     data,\n\u001b[1;32m   1138\u001b[0m     expand_altrep\u001b[38;5;241m=\u001b[39mexpand_altrep,\n\u001b[1;32m   1139\u001b[0m     altrep_constructor_dict\u001b[38;5;241m=\u001b[39maltrep_constructor_dict,\n\u001b[1;32m   1140\u001b[0m     extension\u001b[38;5;241m=\u001b[39mextension,\n\u001b[1;32m   1141\u001b[0m )\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.9/pathlib.py:1259\u001b[0m, in \u001b[0;36mPath.read_bytes\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m   1255\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mread_bytes\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[1;32m   1256\u001b[0m     \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1257\u001b[0m \u001b[38;5;124;03m    Open the file in bytes mode, read it, and close the file.\u001b[39;00m\n\u001b[1;32m   1258\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1259\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[43mmode\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43mrb\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mas\u001b[39;00m f:\n\u001b[1;32m   1260\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m f\u001b[38;5;241m.\u001b[39mread()\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.9/pathlib.py:1252\u001b[0m, in \u001b[0;36mPath.open\u001b[0;34m(self, mode, buffering, encoding, errors, newline)\u001b[0m\n\u001b[1;32m   1246\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mopen\u001b[39m(\u001b[38;5;28mself\u001b[39m, mode\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mr\u001b[39m\u001b[38;5;124m'\u001b[39m, buffering\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m, encoding\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[1;32m   1247\u001b[0m          errors\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, newline\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m   1248\u001b[0m     \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m   1249\u001b[0m \u001b[38;5;124;03m    Open the file pointed by this path and return a file object, as\u001b[39;00m\n\u001b[1;32m   1250\u001b[0m \u001b[38;5;124;03m    the built-in open() function does.\u001b[39;00m\n\u001b[1;32m   1251\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1252\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mio\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mbuffering\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mencoding\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43merrors\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mnewline\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m   1253\u001b[0m \u001b[43m                   \u001b[49m\u001b[43mopener\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_opener\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/anaconda3/lib/python3.9/pathlib.py:1120\u001b[0m, in \u001b[0;36mPath._opener\u001b[0;34m(self, name, flags, mode)\u001b[0m\n\u001b[1;32m   1118\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_opener\u001b[39m(\u001b[38;5;28mself\u001b[39m, name, flags, mode\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0o666\u001b[39m):\n\u001b[1;32m   1119\u001b[0m     \u001b[38;5;66;03m# A stub for the opener argument to built-in open()\u001b[39;00m\n\u001b[0;32m-> 1120\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_accessor\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mopen\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mflags\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmode\u001b[49m\u001b[43m)\u001b[49m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'v11acs_res/v11acs_mag_pos.rda'"
     ]
    }
   ],
   "source": [
    "nm_pos = rdata.parser.parse_file(\"v11acs_nm_pos.rda\")\n",
    "col_pos = rdata.parser.parse_file(\"v11acs_col_pos.rda\")\n",
    "mag_pos = rdata.parser.parse_file(\"v11acs_res/v11acs_mag_pos.rda\")\n",
    "col_err_pos = rdata.parser.parse_file(\"v11acs_res/v11acs_col_err_pos.rda\")\n",
    "\n",
    "nm_pos = rdata.conversion.convert(nm_pos)\n",
    "col_pos = rdata.conversion.convert(col_pos)\n",
    "mag_pos = rdata.conversion.convert(mag_pos)\n",
    "col_err_pos = rdata.conversion.convert(col_err_pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0f06356",
   "metadata": {},
   "outputs": [],
   "source": [
    "nm_pos_list = nm_pos['nm_pos']\n",
    "col_pos_list = col_pos['col_pos']\n",
    "mag_pos_list = mag_pos['mag_pos']\n",
    "col_err_pos_list = col_err_pos['col_err_pos']\n",
    "\n",
    "nm_pos_df = pd.concat(nm_pos_list)\n",
    "col_pos_df = pd.concat(col_pos_list)\n",
    "mag_pos_df = pd.concat(mag_pos_list)\n",
    "col_err_pos_df = pd.concat(col_err_pos_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "id": "2f97d8da",
   "metadata": {},
   "outputs": [],
   "source": [
    "nm_pos_df = nm_pos_df[(nm_pos_df['x'] < 1) & (nm_pos_df['x'] > 0) & (nm_pos_df['y'] < 1) & (nm_pos_df['y'] > 0)]\n",
    "mag_pos_df = mag_pos_df[(mag_pos_df['x'] < 1) & (mag_pos_df['x'] > 0) & (mag_pos_df['y'] < 1) & (mag_pos_df['y'] > 0)]\n",
    "col_pos_df = col_pos_df[(col_pos_df['x'] < 1) & (col_pos_df['x'] > 0) & (col_pos_df['y'] < 1) & (col_pos_df['y'] > 0)]\n",
    "col_err_pos_df = col_err_pos_df[(col_err_pos_df['x'] < 1) & (col_err_pos_df['x'] > 0) & (col_err_pos_df['y'] < 1) & (col_err_pos_df['y'] > 0)]\n",
    "\n",
    "x_nm = np.array(nm_pos_df['x'])\n",
    "y_nm = np.array(nm_pos_df['y'])\n",
    "\n",
    "x_col = np.array(col_pos_df['x'])\n",
    "y_col = np.array(col_pos_df['y'])\n",
    "\n",
    "x_mag = np.array(mag_pos_df['x'])\n",
    "y_mag = np.array(mag_pos_df['y'])\n",
    "\n",
    "x_col_err = np.array(col_err_pos_df['x'])\n",
    "y_col_err = np.array(col_err_pos_df['y'])\n",
    "\n",
    "PDF_nm = fastkde.pdf(x_nm, y_nm, var_names = ['x', 'y'], num_points = 1025)\n",
    "PDF_col = fastkde.pdf(x_col, y_col, var_names = ['x', 'y'], num_points = 1025)\n",
    "PDF_mag = fastkde.pdf(x_mag, y_mag, var_names = ['x', 'y'], num_points = 1025)\n",
    "PDF_col_err = fastkde.pdf(x_col_err, y_col_err, var_names = ['x', 'y'], num_points = 1025)\n",
    "\n",
    "sigma_nm = sum((PDF_nm.values - 1)**2)/512/512\n",
    "sigma_col = sum((PDF_col.values - 1)**2)/512/512\n",
    "sigma_mag = sum((PDF_mag.values - 1)**2)/512/512\n",
    "sigma_col_err = sum((PDF_col_err.values - 1)**2)/512/512"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "cbcb659e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAn8AAAJcCAYAAACIdsJjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAwjklEQVR4nO3df7BldXnn+/fHBkREREUNdLemx2qjjFdS2oI3ZW5QY+gmmtYZqwJYEimdvtyI8VbuqMzEGKNJbow3c41XsKvHoSjyw57EMKZ1WnviWMbJIEk3uYA0DHqCgT42CRdwjMEIdvdz/9gbZ7PZffbq03udfdbe71fVKs5a63vWfs4peHjOs75rfVNVSJIkaT48YdoBSJIkaeVY/EmSJM0Riz9JkqQ5YvEnSZI0Ryz+JEmS5ojFnyRJ0hyx+NPMSfLDSSrJCdOORdLqlmR/kvMneL2/SfKTk7qe1AaLPzXWT2p/l+TJA8feluRLx3G9R5KcMXT85n7x9sPHF7Gk1W7aeaCq/mlVfan/me9P8nttfp60Glj86VidALxzgtf7BnDxoztJ/ifgScu9mN0+qZMmmgckLc3iT8fqw8C/THL6qJNJfizJ3iTf7v/zx8Zc73eBSwf2fw64buiaP53k/03y90kOJHn/wLlHb/G+Nck9wBdHxPTP+92FFzX7ESWtsCXzwFI5oH/+0iR3J3kgyS8P3nrtd/P+MMl1Sb7Tv827aeB7/ybJTybZDPxr4GeT/EOSWwbPD4x/THcwyZsHPvuXhuJ6QpIrk/x1//wfJnn6JH5h0vGw+NOx2gd8CfiXwyf6Se0/Ah8FngH8G+A/JnnGEte7ETgtyQuTrAF+Fhi+7fIQvf8xnA78NPC/JXn90JifAF4IXDAU02XAh4CfrKrbxv94kqZgXB44ag5IcjZwNfAm4EzgqcDaoev/DLCz//27gI8NB1BVnwd+A/j3VXVqVZ0zLuj+Z38ceDNwFr28t25gyC8Ar6eXn84CvgVcNe66Utss/rQc7wPekeSZQ8d/Gvh6Vf1uVR2qqk8C/w143ZjrPfpX/2v64785eLKqvlRVX62qI1V1K/BJesl00Pur6qGq+seBY/878C7g/KpaOIafT9LKO2oeGJMD3gh8pqr+vKoeoZefhhet//Oq2l1Vh/ufM7awa+iNwGer6stV9TDwy8CRgfP/K/BLVbXYP/9+4I1OT9G0+S+gjllV3Zbks8CVwB0Dp84C7h4afjeP/yt82O8CXwY2MHTLFyDJecBvAi8CTgKeCPzR0LADI677LuADVbU45vMlTd9R88CYHHAWA//9V9V3kzwwdO2/Hfj6u8DJSU6oqkPHGfPwZz809NnPBf5DksGC8DDwbIb+yJVWkp0/LdevAP+CxxZ2B+klu0HPYUySq6q76U34vhC4fsSQP6B3q2Z9VT0V2A5k+DIjvu+ngPcm+edLfb6k6RuTB5bKAfcycKs1yZPo3X5dVhgjjj0EnDKw/0MDX98LrB/47FOGPvsAsKWqTh/YTq4qCz9NlcWflqV/G/Xf05vT8qjdwPOTXJLkhCQ/C5wNfLbBJd8KvKqqHhpx7inAg1X1vSTnApc0DHM/sBm4KsnPNPweSdNztDywVA74FPC6/sNmJwG/yuP/OGzq74AfTjL4/8abgYuSnNh/UOSNQ5/92iSv6H/2B3js/1e3A7+e5LkASZ6ZZOsyY5MmxuJPx+MDwA/e+VdVDwCvBf4P4AHg3cBrq+r+cReqqr+uqn1HOf3zwAeSfIfefJ4/bBpgVd3Sj+nfJtnS9Pskrbwl8sBRc0BV7QfeQe+BjnuB7wD3AQ8vI4RHbyU/kOSv+l//MvA8eg9r/Cq9LuTgZ7+9f+ze/pjBaSa/Q69j+Z/6sd8InLeMuKSJStWoLrckSd2T5FTgvwMbq+obUw5HWpXs/EmSOi3J65Kc0l996P8Cvgr8zXSjklav1oq/JNckuS/JyHerpeejSRaS3JrkJW3FIkmaaVvpPXB2ENgIXFTe1pKOqrXbvkn+F+AfgOuq6nErKyS5kN48jQvpzYH4napyLoQkSVKLWuv8VdWXgQeXGLKVXmFYVXUjcHqSM9uKR5IkSdN9yfNaHvti3sX+sXuHBybZBmwDePKTn/zSF7zgBSsSoKTjd9NNN91fVcOrwaxqgzlnDWteegqnTTkiSU19h291LuestGkWf6PewzTyHnRV7QB2AJyWp9fT/mpDm3FJmqibhld9WfWGc855efWUI5LU1BfqU53LOSttmk/7LjLwZnR6b2g/OKVYJEmS5sI0i79dwKX9p35fDny7qh53y1eSJEmT09pt3ySfBM4HzkiySG8t2BMBqmo7vaXALgQW6C20fVlbsUiSJKmnteKvqi4ec77oLYsjSZKkFeIKH5IkSXPE4k+SJGmOWPxJkiTNEYs/SZKkOWLxJ0mSNEcs/iRJkuaIxZ8kSdIcsfiTJEmaIxZ/kiRJc8TiT5IkaY5Y/EmSJM0Riz9JkqQ5YvEnSZI0Ryz+JEmSVqkk1yS5L8ltRzmfJB9NspDk1iQvGXdNiz9JkqTV61pg8xLntwAb+9s24OPjLmjxJ0mStEpV1ZeBB5cYshW4rnpuBE5PcuZS1zxhkgFKkiTNogte+eR64MHDE73mTbc+vB/43sChHVW14xgvsxY4MLC/2D9279G+weJPkiRpjAcePMxf7nnORK+55syvf6+qNh3nZTLiWC31DRZ/kiRJYxRwhCPTDmOURWD9wP464OBS3+CcP0mSpO7aBVzaf+r35cC3q+qot3zBzp8kSVIDxeFa+c5fkk8C5wNnJFkEfgU4EaCqtgO7gQuBBeC7wGXjrmnxJ0mStEpV1cVjzhfw9mO5psWfJEnSGL05f0s+R9EZzvmTJEmaI3b+JEmSGlilT/seMzt/kiRJc8TOnyRJ0hhFcbic8ydJkqSOsfMnSZLUgE/7SpIkqXPs/EmSJI1RwGE7f5IkSeoaO3+SJEkNOOdPkiRJnWPnT5IkaYyCmXnPn8WfJElSA7OxuJu3fSVJkuaKnT9JkqQxivJVL5IkSeoeO3+SJEnjFByejcafnT9JkqR50mrxl2RzkjuTLCS5csT5pyb5TJJbkuxPclmb8UiSJC1H0Xvad5LbtLRW/CVZA1wFbAHOBi5OcvbQsLcDt1fVOcD5wG8nOamtmCRJkuZdm3P+zgUWquougCQ7ga3A7QNjCnhKkgCnAg8Ch1qMSZIkaRnCYTLtICaizdu+a4EDA/uL/WODPga8EDgIfBV4Z1U9rhOaZFuSfUn2fZ+H24pXkgBzjqTZ1mbxN6o8Hn5O5gLgZuAs4EeBjyU57XHfVLWjqjZV1aYTeeKk45SkxzDnSBpWwJGa7DYtbRZ/i8D6gf119Dp8gy4Drq+eBeAbwAtajEmSJGmutVn87QU2JtnQf4jjImDX0Jh7gFcDJHk28CPAXS3GJEmStCyH+/P+JrVNS2sPfFTVoSRXAHuANcA1VbU/yeX989uBDwLXJvkqvdvE76mq+9uKSZIkad61usJHVe0Gdg8d2z7w9UHgp9qMQZIk6XgV+LSvtFrtOXgLew7eMu0wJM0Jc466xuJPM2UwAZuMJbXNnDNfjlQmuk2LxZ9mygVnnTPtECRJWtUs/iRJWib/4Jwfj875m4WnfS3+NHNMxpJWkjlHXdPq077StJiMJa0kc87sK8LhGemZzcZPIUmSpEbs/EmSJDUwzSd0J8nOnyRJ0hyx8ydJkjSGK3xIkiSpk+z8SZIkjRUO12z0zGbjp5AkSVIjdv4kSZLGKODIjPTMZuOnkCRJUiN2/iRJkhrwaV9JkiR1jp0/SZKkMap82leSJEkdZOdPkiSpgSPO+ZMkSVLX2PmTJEkao7e272z0zGbjp5AkSVIjdv4kSZLG8mlfSZIkdZCdP0mSpDFc21eSJEmdZOdPkiSpgcPle/4kSZLUMXb+JEmSxijie/4kSZLUPXb+JEmSGjjie/4kSZLUNXb+JEmSxnBtX0mSJHWSnT9JkqQxivieP0mSJHWPnT9JkqQGXNtXkiRJnWPnT5IkaYwqOOx7/sZLsjnJnUkWklx5lDHnJ7k5yf4kf9ZmPJIkSfOutc5fkjXAVcBrgEVgb5JdVXX7wJjTgauBzVV1T5JntRWPJEnS8oUj+LTvOOcCC1V1V1U9AuwEtg6NuQS4vqruAaiq+1qMR5Ikae61OedvLXBgYH8ROG9ozPOBE5N8CXgK8DtVdd3whZJsA7YBnMwprQQrSY8y50gaVszOnL82i79RvdEa8fkvBV4NPAn4SpIbq+prj/mmqh3ADoDT8vTha0jSRJlzJM2yNou/RWD9wP464OCIMfdX1UPAQ0m+DJwDfA1JkqRVxLV9x9sLbEyyIclJwEXArqExfwL8eJITkpxC77bwHS3GJEmSNNda6/xV1aEkVwB7gDXANVW1P8nl/fPbq+qOJJ8HbgWOAJ+oqtvaikmSJGk5inBkRtb2bfUlz1W1G9g9dGz70P6HgQ+3GYckSdLx8ravJEmSOsfl3SRJksYo4MiMvOplNn4KSZIkNWLnT5Ikaaxw2OXdJEmS1DV2/iRJksZwzp8kSZI6yc6fJElSA875kyRJUufY+ZMkSRqjKs75kyRJUvfY+ZMkSWrgsJ0/SZIkdY2dP0mSpDEKOOLTvpIkSeoaiz9JkqSxwuF6wkS3Rp+abE5yZ5KFJFeOOP/UJJ9JckuS/UkuG3dNiz9JkqRVKMka4CpgC3A2cHGSs4eGvR24varOAc4HfjvJSUtd1zl/kiRJY/TW9l3xOX/nAgtVdRdAkp3AVuD2odCekiTAqcCDwKGlLjq285fkiiRPW27UkrRamM8krTJnJNk3sG0bOr8WODCwv9g/NuhjwAuBg8BXgXdW1ZGlPrRJ5++HgL1J/gq4BthTVdXg+yRptTGfSVq2w5OfLXd/VW1a4vyoVuNwzroAuBl4FfA84E+T/Jeq+vujXXTsT1FV7wU2Av8OeAvw9SS/keR5475XklYT85mkjlkE1g/sr6PX4Rt0GXB99SwA3wBesNRFG5Ww/b+M/7a/HQKeBnwqyW81i12SVgfzmaTlKMKRmuzWwF5gY5IN/Yc4LgJ2DY25B3g1QJJnAz8C3LXURcfe9k3yC8DPAfcDnwDeVVXfT/IE4OvAu5tEL0nTZj6T1CVVdSjJFcAeYA1wTVXtT3J5//x24IPAtUm+Su828Xuq6v6lrttkzt8ZwD+rqruHAjqS5LXL+FkkaVrMZ5KW7cgU3pBXVbuB3UPHtg98fRD4qWO55tjir6ret8S5O47lwyRpmsxnkuR7/iRJksaqgsMr/56/VrjChyRJ0hyx8ydJktTAFFb4aIWdP0mSpDli50+SJGmM3nv+ZqNnNhs/hSRJkhqx8ydJktTA4ZFL7XaPnT9JkqQ5YudPkiRpjMKnfSVJktRBdv4kSZLG8mlfSZIkdZCdP0mSpAaO+LSvJEmSusbOnyRJ0hhVcNinfcdLsjnJnUkWkly5xLiXJTmc5I1txiNJkjTvWuv8JVkDXAW8BlgE9ibZVVW3jxj3IWBPW7FIkiQdL5/2He9cYKGq7qqqR4CdwNYR494B/DFwX4uxSJIkiXbn/K0FDgzsLwLnDQ5IshZ4A/Aq4GVHu1CSbcA2gJM5ZeKBStIgc46kYUVc4aOBUb+hGtr/CPCeqjq81IWqakdVbaqqTSfyxEnFJ0kjmXMkzbI2O3+LwPqB/XXAwaExm4CdSQDOAC5McqiqPt1iXJIkScdsVt7z12bxtxfYmGQD8E3gIuCSwQFVteHRr5NcC3zWwk+SJKk9rRV/VXUoyRX0nuJdA1xTVfuTXN4/v72tz5YkSZqkgpmZ89fqS56rajewe+jYyKKvqt7SZiySJElyhQ9JkqRGfM+fJEmSOsfOnyRJ0jjle/4kSZLUQXb+JEmaoD0Hb3ncsQvOOmcKkWiSitl5z5+dP0mSWrbn4C0ji0JpGiz+JEmaoKN1+ez+dd+R/ry/SW3T4m1fSZImzEJPq5nFnyRJ0hiztMKHt30lSZLmiJ0/SZKkBuz8SZIkqXPs/EmSJI1RzM4KHxZ/kiRJDfiSZ0mSJHWOnT9JkqRxygc+JEmS1EF2/iRJksbwJc+SJEnqJDt/kiRJDdj5kyRJUufY+ZMkSRpjll7ybOdPkiRpjtj5kyRJaqDs/EmSJKlr7PxJkiQ14Nq+kiRJ6hw7f5IkSWOUa/tKkiSpi+z8SZIkNeDTvpIkSeocO3+SJEljucKHJEmSOsjOnyRJUgPO+ZMkSVLn2PmTJEkao/A9f5IkSeogO3+SJEnjVG+Vj1lg50+SJGmO2PmTJElq4AjO+ZMkSVLHtFr8Jdmc5M4kC0muHHH+TUlu7W83JDmnzXgkSZKWo+i952+S27S0VvwlWQNcBWwBzgYuTnL20LBvAD9RVS8GPgjsaCseSZIktTvn71xgoaruAkiyE9gK3P7ogKq6YWD8jcC6FuORJElaJtf2bWItcGBgf7F/7GjeCnxu1Ikk25LsS7Lv+zw8wRAl6fHMOZJmWZudv1Hl8cg35CR5Jb3i7xWjzlfVDvq3hE/L02fkLTuSVitzjqRRZuU9f20Wf4vA+oH9dcDB4UFJXgx8AthSVQ+0GI8kSdLca7P42wtsTLIB+CZwEXDJ4IAkzwGuB95cVV9rMRZJkqTjMs0ndCepteKvqg4luQLYA6wBrqmq/Uku75/fDrwPeAZwdRKAQ1W1qa2YJEmS5l2rK3xU1W5g99Cx7QNfvw14W5sxSJIkHa+q2en8ucKHJEnSHHFtX0mSpAZ8z58kSZI6x86fJElSA7Pynj87f5IkSXPEzp8kSVIDPu0rSZKkzrHzJ0mSNEYRO3+SJEnqHjt/kiRJDczIw752/iRJkuaJnT9JkqRxXNtXkiRJXWTnT5IkqYkZmfRn50+SJGmO2PmTJElqwDl/kiRJalWSzUnuTLKQ5MqjjDk/yc1J9if5s3HXtPMnSZLUQK3wnL8ka4CrgNcAi8DeJLuq6vaBMacDVwObq+qeJM8ad107f5IkSavTucBCVd1VVY8AO4GtQ2MuAa6vqnsAquq+cRe1+JMkSRqj6M35m+QGnJFk38C2behj1wIHBvYX+8cGPR94WpIvJbkpyaXjfhZv+0qSJE3H/VW1aYnzo54wGb75fALwUuDVwJOAryS5saq+drSLWvxJkiSNU8DKP+27CKwf2F8HHBwx5v6qegh4KMmXgXOAoxZ/3vaVJElanfYCG5NsSHIScBGwa2jMnwA/nuSEJKcA5wF3LHVRO3+SJEkNrPTTvlV1KMkVwB5gDXBNVe1Pcnn//PaquiPJ54FbgSPAJ6rqtqWua/EnSZK0SlXVbmD30LHtQ/sfBj7c9JoWf5IkSU3MyNq+Fn+SJEljxeXdJEmS1D12/iRJkpqYkdu+dv4kSZLmiJ0/SZKkcQrn/EmSJKl77PxJkiQ14Zw/SZIkdY2dP0mSpEac8ydJkqSOsfMnSZLUhHP+JEmS1DV2/iRJkpqw8ydJkqSuabX4S7I5yZ1JFpJcOeJ8kny0f/7WJC9pMx5JkqRlKaAy2W1KWiv+kqwBrgK2AGcDFyc5e2jYFmBjf9sGfLyteCRJktRu5+9cYKGq7qqqR4CdwNahMVuB66rnRuD0JGe2GJMkSdKyVE12m5Y2H/hYCxwY2F8EzmswZi1w7+CgJNvodQYBHv5Cfeq2yYa6os4A7p92EMvU5dih2/F3OfbnTjuAY2XOWTW6HDt0O/4ux965nLPS2iz+Rt3MHq5zm4yhqnYAOwCS7KuqTccf3nR0Of4uxw7djr/LsXeROWd16HLs0O34uxx7q3zad6xFYP3A/jrg4DLGSJIkaULaLP72AhuTbEhyEnARsGtozC7g0v5Tvy8Hvl1V9w5fSJIkaepm5Gnf1m77VtWhJFcAe4A1wDVVtT/J5f3z24HdwIXAAvBd4LIGl97RUsgrpcvxdzl26Hb8XY6967r+u+9y/F2OHbodf5dj1xipaT5uIkmS1AFP/OF1deYvvXOi17x727tvmsbcSlf4kCRJmiMWf5qqJC/rr+5ycpInJ9mf5EXTjkvSbDLnaNmqhW1K2nzVizRWVe1Nsgv4NeBJwO9VVZffqSZpFTPnSBZ/Wh0+QO/p8O8BvzDlWCTNPnOOlmG6T+hOkrd9tRo8HTgVeApw8pRjkTT7zDmaaxZ/Wg12AL8M/D7woSnHImn2mXO0PM75k45fkkuBQ1X1B0nWADckeVVVfXHasUmaPeYcyeJPU1ZV1wHX9b8+DJw33YgkzTJzjo7LjLwa2du+kiRJc8TOnyRJUhN2/iRJktQ1dv4kSZLGKXzPnyRJkrrHzp8kSVIDcc6fJEmSusbOnyRJUhN2/iRJktQ1Fn+SJElzxOJPkiRpjjjnT5IkqQGf9pUkSVLn2PmTJElqwhU+JEmS1DV2/iRJksYpfM+fJEmSusfOnyRJUhN2/iRJktQ1dv4kSZIa8D1/kiRJ6hyLP60aSc5PsjjtOCStDuYErTo14W1KLP7UiiSXJNmX5B+S3Jvkc0leMe24JK28rueDJH+T5B/78T+6fWzacUnLZfGniUvyi8BHgN8Ang08B7ga2NriZ65p69qSlm8a+aD/uUvmhCRvSXLtMVzydVV16sB2xVGu+7i59Mean8xnq5idP+nxkjwV+ADw9qq6vqoeqqrvV9VnqupdSZ6Y5CNJDva3jyR54lGu9cIkX0ry35PsT/IzA+euTfLxJLuTPAS8coV+REkNjcsH/TGdzgn9IvK/Jvm/kzwIvH9ULKsxds0viz9N2v8MnAz8h6Oc/yXg5cCPAucA5wLvHR6U5ETgM8B/Ap4FvAP4/SQ/MjDsEuDXgacAfz6Z8CVN0Lh8ALORE84D7urH9esjYvkLVm/saig1+W1aLP40ac8A7q+qQ0c5/ybgA1V1X1X9f8CvAm8eMe7lwKnAb1bVI1X1ReCzwMUDY/6kqv5rVR2pqu9N8GeQNBnj8gF0Jyd8ut+1e3T7FwPnDlbV/1NVh6rqH4djoVfYms+0avieP03aA8AZSU44SsI/C7h7YP/u/rFR4w70E+fg2LUD+weON1hJrRqXD2AFc0KSq+l12ABOAk5I8vr+/j1V9eIlvv31VfWFo5wb9bmDx8xns6Iy7Qgmws6fJu0rwPeA1x/l/EHguQP7z+kfGzVufZInDI395sD+jLxuU5pZ4/IBrGBOqKqfr6rTq+p04OeBP3h0f0zhN86ozx08Zj7TqmLxp4mqqm8D7wOuSvL6JKckOTHJliS/BXwSeG+SZyY5oz/290Zc6i+Ah4B397//fOB1wM4V+UEkHbcG+QDmIyd0OXYNmpGnfb3tq4mrqn+T5O/oTdr+feA7wE30JjP/FXAacGt/+B8BvzbiGo/0n4a7GvhX9P5CvrSq/lv7P4GkSRmTD6D3338XcsJnkhwe2P/TqnpDk29cBbFLj5EqO82SJElLOXnd+lp/xS9O9JoL/+oXb6qqTRO9aAN2/iRJkpqYkX5Za3P+klyT5L4ktx3lfJJ8NMlCkluTvKStWCRJktTT5gMf1wKblzi/BdjY37YBH28xFkmSpOXzJc/jVdWXgQeXGLIVuK56bgROT3JmW/FIkiRpunP+1vLYl1ou9o/dOzwwyTZ63UHWsOalp3DaigQo6fh9h2/dX1XPnHYcx8KcI3VXqzlnRub8TbP4G/Wa7JG/1qraAewAOC1Pr/Py6jbjkjRBX6hP3T1+1OpizpG6q4s5Z6VNs/hbBNYP7K9j9FvdJUmSpm9GOn/TXOFjF3Bp/6nflwPfrqrH3fKVJEnS5LTW+UvySeB8eot6LwK/ApwIUFXbgd3AhcAC8F3gsrZikSRJOl7TfEJ3klor/qrq4jHnC3h7W58vSZKkx5vmbV9JkiStMIs/SZKkOeLavpIkSU3MyJw/O3+SJElzxM6fJEnSOFNej3eS7PxJkiTNETt/kiRJTdj5kyRJUtfY+ZMkSWrCzp8kSZK6xs6fJEnSGMGnfSVJktRBdv4kSZKasPMnSZKkrrHzJ0mSNI4rfEiSJKmL7PxJkiQ1YedPkiRJXWPxJ0mS1ERNeGsgyeYkdyZZSHLlEuNeluRwkjeOu6bFnyRJ0iqUZA1wFbAFOBu4OMnZRxn3IWBPk+ta/EmSJDWQmuzWwLnAQlXdVVWPADuBrSPGvQP4Y+C+Jhe1+JMkSZqOM5LsG9i2DZ1fCxwY2F/sH/uBJGuBNwDbm36oT/tKkiQ1Mfmnfe+vqk1LnE+DKD4CvKeqDiejhj+exZ8kSdLqtAisH9hfBxwcGrMJ2Nkv/M4ALkxyqKo+fbSLWvxJkiSNcwxP6E7QXmBjkg3AN4GLgEseE1bVhke/TnIt8NmlCj+w+JMkSVqVqupQkivoPcW7BrimqvYnubx/vvE8v0EWf5IkSQ1MY23fqtoN7B46NrLoq6q3NLmmT/tKkiTNETt/kiRJTbi2ryRJkrrGzp8kSVID05jz1wY7f5IkSXPEzp8kSVITdv4kSZLUNXb+JEmSxpnOCh+tsPMnSZI0R+z8SZIkjZH+Ngvs/EmSJM0RO3+SJElNOOdPkiRJXdNq8Zdkc5I7kywkuXLE+acm+UySW5LsT3JZm/FIkiQtV2qy27S0VvwlWQNcBWwBzgYuTnL20LC3A7dX1TnA+cBvJzmprZgkSZLmXZudv3OBhaq6q6oeAXYCW4fGFPCUJAFOBR4EDrUYkyRJ0vLUhLcpabP4WwscGNhf7B8b9DHghcBB4KvAO6vqyPCFkmxLsi/Jvu/zcFvxShJgzpE029os/ka9Dme4zr0AuBk4C/hR4GNJTnvcN1XtqKpNVbXpRJ446Tgl6THMOZJGsvM31iKwfmB/Hb0O36DLgOurZwH4BvCCFmOSJEmaa20Wf3uBjUk29B/iuAjYNTTmHuDVAEmeDfwIcFeLMUmSJB27CT/pO82nfVt7yXNVHUpyBbAHWANcU1X7k1zeP78d+CBwbZKv0rtN/J6qur+tmCRJkuZdqyt8VNVuYPfQse0DXx8EfqrNGCRJkibCFT4kSZLUNa7tq5mz5+AtAFxw1jlTjkSSNEumOU9vkuz8aWY9WgRKUpv2HLzlB5vUBRZ/mmkmY0nSxPieP2n189avpLaZZ+aHr3qRVikTsaSVZt5Rl1j8SZIkjTPlW7WT5G1fSZKkOWLnT5IkqQk7f5IkSeoaO3+SJEljBF/yLEmSpA6y8ydJktSEnT9JkiR1jZ0/SZKkBlKz0fqz8ydJkjRH7PxJkiSN4wofkiRJ6iI7f5IkSQ34nj9JkiR1jp0/SZKkJuz8SZIkqWvs/EmSJDXgnD9JkiR1jp0/SZKkJuz8SZIkqWvs/EmSJI1TzvmTJElSB9n5kyRJasLOnyRJkrrGzp8kSdIYwTl/kiRJ6iA7f5IkSU3UbLT+7PxJkiTNETt/kiRJDTjnT5IkSZ1j50+SJGmcwvf8SZIkqXvs/EmSJDWQI9OOYDLs/EmSJM0RO3+SJElNOOdvvCSbk9yZZCHJlUcZc36Sm5PsT/JnbcYjSZI071rr/CVZA1wFvAZYBPYm2VVVtw+MOR24GthcVfckeVZb8UiSJB0P3/M33rnAQlXdVVWPADuBrUNjLgGur6p7AKrqvhbjkSRJmnttFn9rgQMD+4v9Y4OeDzwtyZeS3JTk0lEXSrItyb4k+77Pwy2FK0k95hxJj1P01vad5DYlbT7wkRHHhn/SE4CXAq8GngR8JcmNVfW1x3xT1Q5gB8BpefqMNF0lrVbmHEmzrM3ibxFYP7C/Djg4Ysz9VfUQ8FCSLwPnAF9DkiRpFXHO33h7gY1JNiQ5CbgI2DU05k+AH09yQpJTgPOAO1qMSZIkaa611vmrqkNJrgD2AGuAa6pqf5LL++e3V9UdST4P3AocAT5RVbe1FZMkSdKyzUjnr9WXPFfVbmD30LHtQ/sfBj7cZhySJEnqcYUPSZKkMYJz/iRJktRBdv4kSZLGmfK7+SbJzp8kSdIcsfMnSZLUgHP+JEmS1Dl2/iRJkpqw8ydJkqSusfMnSZLUgHP+JEmS1Dl2/iRJksYp4MhstP7s/EmSJM0RO3+SJElNzEbjz86fJEnSPLHzJ0mS1IBP+0qSJKlz7PxJkiQ1UbPR+hvb+UtyRZKnrUQwktQm85mk45Ga7DYtTW77/hCwN8kfJtmcJG0HJUktMZ9Jmntji7+qei+wEfh3wFuAryf5jSTPazk2SZoo85mkZasWtilp9MBHVRXwt/3tEPA04FNJfqvF2CRp4sxnkrqkf5fiziQLSa4ccf5NSW7tbzckOWfcNcc+8JHkF4CfA+4HPgG8q6q+n+QJwNeBdx/7jyJJK898Jmm5AmSFH/hIsga4CngNsEhv2squqrp9YNg3gJ+oqm8l2QLsAM5b6rpNnvY9A/hnVXX34MGqOpLktcfyQ0jSlJnPJHXJucBCVd0FkGQnsBX4QfFXVTcMjL8RWDfuomOLv6p63xLn7hj3/ZK0WpjPJB2XIxO/4hlJ9g3s76iqHQP7a4EDA/uLLN3VeyvwuXEf6nv+JEmSpuP+qtq0xPlRbyQYee85ySvpFX+vGPehFn+SJEkNrPScP3qdvvUD++uAg8ODkryY3jzmLVX1wLiLurybJEnS6rQX2JhkQ5KTgIuAXYMDkjwHuB54c1V9rclF7fxJkiSNM4V381XVoSRXAHuANcA1VbU/yeX989uB9wHPAK7uv7f+0JhbyRZ/kiRJq1VV7QZ2Dx3bPvD124C3Hcs1Lf4kSZLGKlj5OX+tcM6fJEnSHLHzJ0mS1EBmo/Fn50+SJGme2PmTJElqwjl/kiRJ6ho7f5IkSeMUZPJr+06FnT9JkqQ5YudPkiSpCef8SZIkqWvs/EmSJDUxG42/djt/STYnuTPJQpIrlxj3siSHk7yxzXgkSZLmXWudvyRrgKuA1wCLwN4ku6rq9hHjPgTsaSsWSZKk4xXn/I11LrBQVXdV1SPATmDriHHvAP4YuK/FWCRJkkS7xd9a4MDA/mL/2A8kWQu8Adi+1IWSbEuyL8m+7/PwxAOVpEHmHEkjVU12m5I2i7+MODb8k34EeE9VHV7qQlW1o6o2VdWmE3nipOKTpJHMOZJmWZtP+y4C6wf21wEHh8ZsAnYmATgDuDDJoar6dItxSZIkHZsCZmSFjzaLv73AxiQbgG8CFwGXDA6oqg2Pfp3kWuCzFn6SJEntaa34q6pDSa6g9xTvGuCaqtqf5PL++SXn+UmSJK0WoWbmad9WX/JcVbuB3UPHRhZ9VfWWNmORJEmSK3xIkiQ1MyOdP9f2lSRJmiN2/iRJkpqYkc6fxZ/mwp6Dt4w8fsFZ56xwJJLmyZ6Dt5hntOp421cz72iF37hzkrRcew7e8oP8Yp6ZEY++52+S25RY/EmSJM0Rb/tq5g3echn+C9zbMZLaYG6ZTb7nT+ogE7Ikad5Z/EmSJDUxI50/5/xJkiTNETt/kiRJY5WdP0mSJHWPnT9JkqRxCjt/kiRJ6h47f5IkSU1McVWOSbLzJ0mSNEfs/EmSJDUwKyt82PmTJEmaI3b+JEmSmrDzJ0mSpK6x8ydJkjROAUfs/EmSJKlj7PxJkiSN5dq+kiRJ6iA7f5IkSU3Y+ZMkSVLX2PmTJElqws6fJEmSusbOnyRJ0ji+50+SJEldZOdPkiRprII6Mu0gJsLiT5IkqQkf+JAkSVLX2PmTJEkaxwc+JEmS1EV2/iRJkppwzp8kSZK6xs6fJElSE3b+JEmS1DV2/iRJksYqO39NJNmc5M4kC0muHHH+TUlu7W83JDmnzXgkSZLmXWudvyRrgKuA1wCLwN4ku6rq9oFh3wB+oqq+lWQLsAM4r62YJEmSlqWAI7OxvFubnb9zgYWququqHgF2AlsHB1TVDVX1rf7ujcC6FuORJEmae23O+VsLHBjYX2Tprt5bgc+NOpFkG7AN4GROmVR8kjSSOUfSSDMy56/N4i8jjo38rSV5Jb3i7xWjzlfVDnq3hDktT5+N37ykVcucI2mWtVn8LQLrB/bXAQeHByV5MfAJYEtVPdBiPJIkScs3I52/Nuf87QU2JtmQ5CTgImDX4IAkzwGuB95cVV9rMRZJkiTRYuevqg4luQLYA6wBrqmq/Uku75/fDrwPeAZwdRKAQ1W1qa2YJEmSlqfgyGx0/lp9yXNV7QZ2Dx3bPvD124C3tRmDJEmS/gdX+JAkSRqnoMr3/EmSJKlj7PxJkiQ1MSNz/uz8SZIkzRE7f5IkSU34nj9JkiR1jZ0/SZKkcargiE/7SpIkqWPs/EmSJDXhnD9JkiR1jZ0/SZKkBso5f5IkSeoaO3+SJEljlXP+JEmS1D12/iRJksYpXNtXkiRJ3WPnT5IkqYnyaV9JkiR1jJ0/SZKkMQoo5/xJkiSpa+z8SZIkjVPlnD9JkiR1j50/SZKkBpzzJ0mSpFYl2ZzkziQLSa4ccT5JPto/f2uSl4y7pp0/SZKkJlZ4zl+SNcBVwGuARWBvkl1VdfvAsC3Axv52HvDx/j+Pys6fJEnS6nQusFBVd1XVI8BOYOvQmK3AddVzI3B6kjOXuqidP0mSpDG+w7f2fKE+dcaEL3tykn0D+zuqasfA/lrgwMD+Io/v6o0asxa492gfavEnSZI0RlVtnsLHZsSx4adOmox5DG/7SpIkrU6LwPqB/XXAwWWMeQyLP0mSpNVpL7AxyYYkJwEXAbuGxuwCLu0/9fty4NtVddRbvuBtX0mSpFWpqg4luQLYA6wBrqmq/Uku75/fDuwGLgQWgO8Cl427rsWfJEnSKlVVu+kVeIPHtg98XcDbj+Wa3vaVJEmaIxZ/kiRJc8TiT5IkaY5Y/EmSJM0Riz9JkqQ5YvEnSZI0Ryz+JEmS5ojFnyRJ0hyx+JMkSZojrRZ/STYnuTPJQpIrR5xPko/2z9+a5CVtxiNJkjTvWiv+kqwBrgK2AGcDFyc5e2jYFmBjf9sGfLyteCRJktTu2r7nAgtVdRdAkp3AVuD2gTFbgev669LdmOT0JGdW1b1Hu+jzX/pP+NN9f9Ri2JImKcm0Q5AkDWiz+FsLHBjYXwTOazBmLfCY4i/JNnqdQYCHk9w22VBX1BnA/dMOYpm6HDt0O/4ux/7caQdwrIZzzhfqU+ac6ehy7NDt+Lsce+dyzkprs/gb9ed+LWMMVbUD2AGQZF9VbTr+8Kajy/F3OXbodvxdjr2LzDmrQ5djh27H3+XYNV6bD3wsAusH9tcBB5cxRpIkSRPSZvG3F9iYZEOSk4CLgF1DY3YBl/af+n058O2l5vtJkiTp+LR227eqDiW5AtgDrAGuqar9SS7vn98O7AYuBBaA7wKXNbj0jpZCXildjr/LsUO34+9y7F3X9d99l+PvcuzQ7fi7HLvGSO9BW0mSJM0DV/iQJEmaIxZ/kiRJc2TVFn9dXxquQfxv6sd9a5IbkpwzjThHGRf7wLiXJTmc5I0rGd84TeJPcn6Sm5PsT/JnKx3j0TT49+apST6T5JZ+7E3myaqBLuecLucb6HbO6XK+AXPO3KqqVbfRe0Dkr4F/ApwE3AKcPTTmQuBz9N4V+HLgL6Yd9zHG/2PA0/pfb1kt8TeJfWDcF+k9tPPGacd9jL/70+mtNPOc/v6zph33McT+r4EP9b9+JvAgcNK0Y+/61uWc0+V80zT+gXGrKud0Od8cQ/zmnBncVmvn7wdLw1XVI8CjS8MN+sHScFV1I3B6kjNXOtCjGBt/Vd1QVd/q795I7x2Hq0GT3z3AO4A/Bu5byeAaaBL/JcD1VXUPQFWtlp+hSewFPCW9NdNOpZeID61smDOpyzmny/kGup1zupxvwJwzt1Zr8Xe0Zd+Odcy0HGtsb6XXUVgNxsaeZC3wBmD7CsbVVJPf/fOBpyX5UpKbkly6YtEtrUnsHwNeSO9l6F8F3llVR1YmvJnW5ZzT5XwD3c45Xc43YM6ZW20u73Y8JrY03JQ0ji3JK+kl41e0GlFzTWL/CPCeqjrc+2NwVWkS/wnAS4FXA08CvpLkxqr6WtvBjdEk9guAm4FXAc8D/jTJf6mqv285tlnX5ZzT5XwD3c45Xc43YM6ZW6u1+Ov60nCNYkvyYuATwJaqemCFYhunSeybgJ39JHwGcGGSQ1X16RWJcGlN/925v6oeAh5K8mXgHGDaybhJ7JcBv1lVBSwk+QbwAuAvVybEmdXlnNPlfAPdzjldzjdgzplf0550OGqjV5TeBWzgf0xC/adDY36ax06+/stpx32M8T+H3somPzbteI819qHx17JKJl8fw+/+hcB/7o89BbgNeFFHYv848P7+188GvgmcMe3Yu751Oed0Od80jX9o/KrJOV3ON8cQvzlnBrdV2fmr9paGWxEN438f8Azg6v5fs4eqatO0Yn5Uw9hXrSbxV9UdST4P3AocAT5RVbdNL+qehr/7DwLXJvkqvSLkPVV1/9SCnhFdzjldzjfQ7ZzT5XwD5px55vJukiRJc2S1Pu0rSZKkFlj8SZIkzRGLP0mSpDli8SdJkjRHLP4kSZLmiMWfJEnSHLH4kyRJmiMWf5qqJC9LcmuSk5M8Ocn+JC+adlySZpM5R/Ilz1oFkvwacDK9Rc8Xq+r/nHJIkmaYOUfzzuJPU5fkJGAv8D16a48ennJIkmaYOUfzztu+Wg2eDpwKPIXeX+OS1CZzjuaanT9NXZJdwE5gA3BmVV0x5ZAkzTBzjubdCdMOQPMtyaXAoar6gyRrgBuSvKqqvjjt2CTNHnOOZOdPkiRprjjnT5IkaY5Y/EmSJM0Riz9JkqQ5YvEnSZI0Ryz+JEmS5ojFnyRJ0hyx+JMkSZoj/z9jqiluCOpyfwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "PDF_nm.values = PDF_nm.values > 5*sigma_nm\n",
    "PDF_col.values = PDF_col.values > 0.4\n",
    "PDF_mag.values = PDF_mag.values > 0.4\n",
    "PDF_col_err.values = PDF_col_err.values > 0.4\n",
    "\n",
    "# Assuming PDF_nm, PDF_mag, PDF_col, PDF_col_err are xarray DataArray objects\n",
    "fig, axs = plt.subplots(2, 2, figsize=(10, 10), sharex=True, sharey=True)\n",
    "\n",
    "# Define a normalization for the color scale\n",
    "norm = mcolors.Normalize(vmin=0, vmax=1)\n",
    "\n",
    "# Plot each DataArray with the same normalization\n",
    "im1 = PDF_nm.plot(ax=axs[0, 0], add_colorbar=False, norm=norm)\n",
    "im2 = PDF_mag.plot(ax=axs[0, 1], add_colorbar=False, norm=norm)\n",
    "im3 = PDF_col.plot(ax=axs[1, 0], add_colorbar=False, norm=norm)\n",
    "im4 = PDF_col_err.plot(ax=axs[1, 1], add_colorbar=False, norm=norm)\n",
    "\n",
    "fig.subplots_adjust(right=0.85)\n",
    "cbar_ax = fig.add_axes([0.88, 0.15, 0.03, 0.7])  # [left, bottom, width, height]\n",
    "\n",
    "axs[0,0].set_title('No Mark')\n",
    "axs[0,1].set_title('Magnitude')\n",
    "axs[1,0].set_title('Color')\n",
    "axs[1,1].set_title('Color + Error')\n",
    "\n",
    "\n",
    "# Create a colorbar in the specified axis\n",
    "cbar = fig.colorbar(im1, cax=cbar_ax)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3685b49",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/davidolohowski/anaconda3/lib/python3.9/site-packages/fastkde/fastKDE.py:1154: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.\n",
      "  inputVariables = npy.array(var0[npy.newaxis, :])\n",
      "/Users/davidolohowski/anaconda3/lib/python3.9/site-packages/fastkde/fastKDE.py:1180: FutureWarning: Support for multi-dimensional indexing (e.g. `obj[:, None]`) is deprecated and will be removed in a future version.  Convert to a numpy array before indexing instead.\n",
      "  varn = npy.array(var_args[i][npy.newaxis, :])\n"
     ]
    }
   ],
   "source": [
    "np.random.seed(537380)\n",
    "nm_sim_list = []\n",
    "mag_sim_list = []\n",
    "col_sim_list = []\n",
    "col_err_sim_list = []\n",
    "for j in range(100):\n",
    "    block_idx = np.random.choice(range(1800), 1800, replace=True)\n",
    "\n",
    "    # Create resampled blocks\n",
    "    rs_blocks = [range(100 * x, 100 * (x+1)) for x in block_idx]\n",
    "    rs_blocks_flat = [item for sublist in rs_blocks for item in sublist]\n",
    "\n",
    "    X_nm_j = pd.concat([nm_pos_list[i] for i in rs_blocks_flat], ignore_index=True)\n",
    "    X_nm_j = X_nm_j[(X_nm_j['x'] < 1) & (X_nm_j['x'] > 0) & (X_nm_j['y'] < 1) & (X_nm_j['y'] > 0)]\n",
    "\n",
    "    X_mag_j = pd.concat([mag_pos_list[i] for i in rs_blocks_flat], ignore_index=True)\n",
    "    X_mag_j = X_mag_j[(X_mag_j['x'] < 1) & (X_mag_j['x'] > 0) & (X_mag_j['y'] < 1) & (X_mag_j['y'] > 0)]\n",
    "    \n",
    "    X_col_j = pd.concat([col_pos_list[i] for i in rs_blocks_flat], ignore_index=True)\n",
    "    X_col_j = X_col_j[(X_col_j['x'] < 1) & (X_col_j['x'] > 0) & (X_col_j['y'] < 1) & (X_col_j['y'] > 0)]\n",
    "    \n",
    "    X_col_err_j = pd.concat([col_err_pos_list[i] for i in rs_blocks_flat], ignore_index=True)\n",
    "    X_col_err_j = X_col_err_j[(X_col_err_j['x'] < 1) & (X_col_err_j['x'] > 0) & (X_col_err_j['y'] < 1) & (X_col_err_j['y'] > 0)]\n",
    "\n",
    "\n",
    "    PDF_nm_sim = fastkde.pdf(X_nm_j['x'], X_nm_j['y'], var_names = ['x', 'y'], num_points = 1025)\n",
    "    PDF_mag_sim = fastkde.pdf(X_mag_j['x'], X_mag_j['y'], var_names = ['x', 'y'], num_points = 1025)\n",
    "    PDF_col_sim = fastkde.pdf(X_col_j['x'], X_col_j['y'], var_names = ['x', 'y'], num_points = 1025)\n",
    "    PDF_col_err_sim = fastkde.pdf(X_col_err_j['x'], X_col_err_j['y'], var_names = ['x', 'y'], num_points = 1025)\n",
    "\n",
    "    \n",
    "    nm_pdf = PDF_nm_sim.values.tolist()\n",
    "    nm_sim_list.append(nm_pdf)\n",
    "    mag_pdf = PDF_mag_sim.values.tolist()\n",
    "    mag_sim_list.append(mag_pdf)\n",
    "    col_pdf = PDF_col_sim.values.tolist()\n",
    "    col_sim_list.append(col_pdf)\n",
    "    col_err_pdf = PDF_col_err_sim.values.tolist()\n",
    "    col_err_sim_list.append(col_err_pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "a20c9630",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "range(1, 100)\n"
     ]
    }
   ],
   "source": [
    "print(range(1, 100))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a30c1115",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
